if (!require("tidyverse")) install.packages("tidyverse")
if (!require("rms")) install.packages("rms")
if (!require("rmda")) install.packages("rmda")
if (!require("ggDCA")) devtools::install_github("yikeshu0611/ggDCA")
if (!require("ggplot2")) install.packages("ggplot2")

rm(list = ls())

df_select <- readRDS("~/analysis/lyz_ml/rds/step_03_mlr3_select.RDS")

colnames(df_select) <- make.names(names(df_select))

df_select$group <- as.numeric(df_select$group) - 1

# ggplot 画决策曲线
if (T) {
  ddist <- datadist(df_select)
  options(datadist = "ddist")

  all_factors <- lrm(group ~ ., df_select)
  #  ggDCA::dca(all_factors) %>% ggplot()

  bronchoscopy <- lrm(group ~ bronchoscopy, df_select)
  age.years <- lrm(group ~ age.years, df_select)
  days.of.illness.prior.to.admission <- lrm(group ~ days.of.illness.prior.to.admission, df_select)
  wbc.before.treatment <- lrm(group ~ wbc.before.treatment, df_select)
  plr.before.treatment <- lrm(group ~ plr.before.treatment, df_select)
  crp.before.treatment <- lrm(group ~ crp.before.treatment, df_select)
  ldh.before.treatment <- lrm(group ~ ldh.before.treatment, df_select)

  # 模型合并比较

  # all factor
  Model_data <- dca(all_factors, bronchoscopy,age.years, days.of.illness.prior.to.admission, wbc.before.treatment,
                    plr.before.treatment,crp.before.treatment, ldh.before.treatment,
                    model.names = c("All factors", 'Bronchoscopy', 'Age', 'Days of illness prior to admission', 'WBC before treatment',
                                    'PLR before treatment',"CRP before treatment", 'LDH before treatment')
  )

  rbind(Model_data) %>%
    ggplot() +
    geom_line(aes(color = model)) +
    ggsci::scale_color_d3(palette = "category20") +
    scale_x_continuous( breaks = seq(0.1,1,0.2),limits = c(0.1,0.9))+
    theme(
      strip.text.x = element_blank(),
      strip.text.y = element_blank()
    )

  ggsave(filename = "result/step_07_dca.pdf", width = 32, height = 20, dpi = 300, units = "cm")
}

# 默认决策曲线
if (F) {
  Models <- decision_curve(
    group ~ pda +
      nec + age_days +
      birth_weight + vsd_shunting +
      gestation_age + admission_weight +
      pulmory_hemorrhage + pda_maximum_diameter +
      pda_blood_flow_velocity + pda_blood_flow_direction +
      area_of_tricuspid_regurgitation + peak_velocity_of_tricuspid_regurgitation,
    data = df_select,
    family = binomial(link = "logit"),
    thresholds = seq(0, 1, by = 0.01),
    confidence.intervals = 0.95,
    study.design = "case-control",
    population.prevalence = 0.3
  )

  plot_decision_curve(Models,
    cost.benefit.axis = FALSE,
    col = c("red2", "blue2", "pink2"),
    curve.names = c("Model"),
    confidence.intervals = FALSE,
    standardize = FALSE
  )

  pdf(file = "result/step_07_dca_impact.pdf", width = 12, height = 6)
  plot_clinical_impact(Models,
    population.size = 1000,
    cost.benefit.axis = T,
    n.cost.benefits = 8,
    col = c("red3", "blue3", "pink3"),
    ylim = c(0, 1000),
    xlim = c(0, 1),
    legend.position = "topright",
    # legend.position = "none",
    confidence.intervals = T,
    # cost.benefit.xlab = "Cost:Benefit Ratio",
    cost.benefit.xlab = "Benefit Ratio",
    ylab = "Number high risk",
    xlab = "High Rish Threshold",
    lwd = 3,
    lty = 2
  )
  dev.off()

  # 多个模型展示
  if (F) {
    pda <- decision_curve(
      group ~ pda,
      data = df_select,
      family = binomial(link = "logit"),
      thresholds = seq(0, 1, by = 0.01),
      confidence.intervals = 0.95, study.design = "case-control",
      population.prevalence = 0.3
    )

    List <- list(Models, pda)

    plot_decision_curve(List,
      cost.benefit.axis = FALSE,
      col = c("red", "blue"),
      # "topright", "right", "bottomright", "bottom",
      # "bottomleft", "left", "topleft", "top", "none"
      legend.position = "right",
      xlim = c(0, 1),
      curve.names = c("Model", "PDA"),
      confidence.intervals = FALSE,
      standardize = FALSE
    )
  }
}
